#ifndef _SOLVER_H
#define _SOLVER_H
#include "Geometry.h"
#include "Partition.h"
#include "Excitation.h"

class Solver
{
public:
	Solver(){};
	~Solver(){};
};
class MoM_Solver:public Solver
{
public:
	MoM_Solver(MoM_Partition* part)
		:Geo(part->Geo),Part(part),Ex(part->Ex)
	{
		AllocateMemory();
		FillZ();
		for (int i=0; i<Part->EdgeSet.size(); i++)
		{
			Part->Pre(i) = 1.0/sqrt(Part->Z(i,i));
		}
		for (int i=0; i<Part->EdgeSet.size(); i++)
		{
			for (int j=0; j<Part->EdgeSet.size(); j++) 
			{
				Part->Z(i,j) = Part->Pre(i) * Part->Z(i,j) * Part->Pre(j);
			}
		}
	};
	
	~MoM_Solver(){};
	void Execution(real_data eps)
	{
		FillV();
		Part->V *= Part->Pre;
		CG(Part->EdgeSet.size(), Part->V, Itmax, eps);
// 		GMRES(Part->EdgeSet.size(), 30, Part->V, eps);
		Part->I *= Part->Pre;
		for (int i=0; i<Part->Edge_Num; i++)
		{
			Geo->I(Part->Edge[i]) = Part->I(i);
		}
	};
private:
	void FillZ();

	void AllocateMemory();

	void Rules(int Node1Fld, int Node2Fld, int Node1Src, int Node2Src, bool& NearEdge,
		int& RuleFld, int& RuleSrc);

	void Rules2(int PatFld, int PatSrc, bool& NearPat, int& RuleFld, int& RuleSrc);

	dComplex CalcMtx(const int FldEdge,const int SrcEdge);

	int IGetN3(int Pat, int Node1, int Node2)
	{
		return Geo->Patch(Pat)[0]+Geo->Patch(Pat)[1]+Geo->Patch(Pat)[2]-Node1-Node2;
	}

	dComplex Cfune(real_data RiRj, const State& State);

	dComplex Cfunh(const State& sState);

	real_data Rfunh(const Position3D& FldNorm, const Position3D& RhoFld, const SingularInt& SignInt);

	void Proj(const Position3D& R, const Position3D& R1, const Position3D& R2
		, const Position3D& R3, const Position3D& R30, Position3D& SrcNorm
		, real_data Aera, SingularInt& SingInt,  State& sState);

	void Trans(const Position3D& Bn, const Position3D& R, const Position3D& R0
		,  const Position3D& R2, const Position3D& R3, Position3D& U1
		, real_data& P10, real_data& R10, real_data& R11, blitz::TinyVector<real_data,2>& Al1);

	void FillV();

	void CG(const int Unknown, const blitz::Array<dComplex,1>& Cb, const int Itmax, const real_data eps);

	void GMRES(const int iUnknown, const int M, const blitz::Array<dComplex,1>& Cb, const real_data eps);

private:
	MoM_Geometry* Geo;
	MoM_Partition* Part;
	const PlaneWave* Ex;
};

#endif